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Abstract 

An in situ computer vision method is presented 
for omnidirectional bee counting on landing pads of 
Langstroth beehives used by many apiarists world- 
wide. Bee counts are computed from static images 
of beehive landing pads. The presented method au- 
tomates landing pad localization and improves bee 
counts through elimination of landing pad skew. Two 
bee counting algorithms are presented for counting 
bees on localized landing pads. The first bee counting 
algorithm is based on the ID Haar Wavelet Trans- 
form. The second algorithm is based on contour anal- 
ysis. The relative performance of both algorithms is 
evaluated on a sample of 793 images obtained from 
two electronic beehive monitoring devices deployed in 
live Langstorth beehives in northern Utah. 

Keywords: Wavelets, Haar Wavelet Transform, 
Contour Analysis, Electronic Beehive Monitoring, 
Sustainable Computing 

Nomenclature: EBM - Electronic Beehive Mon- 
itoring, CV - Computer Vision, EBMD - Electronic 
Beehive Monitoring Device, RPi - Raspberry Pi, 
HWT - Haar Wavelet Transform 


1 Introduction 


Professional and amateur beekeepers use visual esti- 
mates of forager traffic levels to evaluate bee colonies’ 
health. Higher levels of forager traffic may indicate 
onsets of swarms or colony robbing activities; lower 
levels of forager traffic may indicate mite infestations, 
failing queens, or chemical poisonings [1]. Continuous 
and rapid advances in electronic sensors have made 
it feasible to transform individual beehives and entire 
apiaries into ad hoc sensor networks that monitor bee 
colonies’ health in situ [2]. Electronic beehive mon- 
itoring (EBM) can also help researchers and practi- 
tioners to collect critical data on colony behavior and 
phenology without invasive beehive inspections [3]. 


In this article, an in situ computer vision (CV) 
method is presented for omnidirectional bee counting 
on landing pads of Langstroth beehives [4] used by 
many apiarists all over the world. Figure (1) shows 
the flow chart of the vision-based bee counting method 
realized in the software system described in this arti- 
cle. 
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Figure 1: Flow chart of the method 


As shown in figure (1), the system takes as input 
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a static image of a landing pad. The image is put 
through the pad localization module to localize the 
landing pad. Since the pad may be skewed, the image 
of the localized pad is processed by the skew detection 
module where the pad’s skew angle is determined and 
the image is rotated accordingly. The localized and 
rotated pad is given to the bee counting module to 
compute a bee count, a non-negative real number ap- 
proximating the number of bees on the pad. 

The presented bee counting method is omnidirec- 
tional in that it does not distinguish between incoming 
and outgoing bee traffic: in situ bee counts are com- 
puted from static images of landing pads regardless 
of the actual orientations of the bees. Nor does the 
algorithm currently distinguish among three different 
types of bees: the worker, the drone, and the queen. 
The distinction of the three bee classes is beyond the 
current hardware capabilities of the presented EBM 
system in that it will likely require higher resolution 
cameras placed closer to landing pads. The identifica- 
tion of the queen bee on the landing pad is not a seri- 
ous limitation, because, under normal circumstances, 
the mated queen never appears on the landing pad. 

A fundamental objective of the reported research 
is to argue by demonstration that, given increasing 
availability and decreasing costs of CV hardware (e.g., 
miniature cameras and storage devices) and a prolif- 
eration of free, reliable open-source CV tools and li- 
braries, CV is ready for inclusion in standard EBM 
sensor suites. 

Another fundamental, long-term objective of this 
project is to create a suite of replicable hardware and 
software tools for citizen scientists to build their own 
EBM devices (EBMDs) and to promote the grassroots 
development of regional, national, and international 
cyberinfrastructures for capturing, sharing, and an- 
alyzing EBM data [5]. The hardware of the EBM 
system described in this article was assembled with 
off-the-shelf hardware components. The system’s soft- 
ware was developed and tested with open source soft- 
ware tools. 

The remainder of this article is organized as fol- 
lows. In Section (2), related work is reviewed. In Sec- 
tion (3), the hardware and software details of BeePi 

[6] , a multi-sensor, low-energy, solar-powered EBMD, 
are presented. Section (4) presents a method for lo- 
calizing landing pads in images captured by BeePi. In 
Section (5), two bee counting algorithms are presented 
to count bees on localized landing pads. Section (6) 
describes an evaluation of both algorithms on a sam- 
ple of 793 images and analyzes the results. In Section 

(7) , conclusions are drawn. 

2 Related Work 

EBM is a mature research and development topic. In 
the 1950’s, Woods, a sound engineer, designed and 
built Apidictor, an audio beehive monitoring tool [7]. 


Bencsik [8] equipped several hives with accelerometers 
and observed increasing amplitudes a few days before 
swarming, with a sharp change at the point of swarm- 
ing. Evans [9] developed Arnia, a smartphone-based 
beehive monitoring system that uses weight, temper- 
ature, humidity, and sound. The system breaks down 
hive sounds into flight buzzing, fanning, and ventilat- 
ing and sends digital alerts to beekeepers. 

Ferrari et al. [10] assembled a system for monitor- 
ing swarm sounds in beehives. The system included 
a microphone, a temperature sensor, and a humidity 
sensor placed in a beehive and connected via under- 
ground cables to a computer in a nearby barn. 

Rangel and Seeley [11] investigated signals of hon- 
eybee swarms. Five custom designed observation 
hives were sealed with glass covers. The captured 
video and audio data were monitored daily by human 
observers. The researchers found that approximately 
one hour before swarm exodus, the production of pip- 
ing signals gradually increased and ultimately peaked 
at the start of the swarm departure. 

In a longitudinal field investigation, Meikle and 
Holst [12] placed four beehives on precision electronic 
scales linked to data loggers to record weight for 
over sixteen months. The researchers investigated 
the effect of swarming on daily data and reported 
that empty beehives also had detectable daily weight 
changes due to moisture level changes in the wood. 

Bromenshenk et al. [13] developed electronic 
SmartHive devices equipped with electronic scales, in- 
frared bee counters, temperature and humidity sen- 
sors, digital weather stations, and wireless communi- 
cation lines for Internet-based remote monitoring. 

Silva et al. [14] proposed an intelligent trap with 
a laser sensor to selectively classify various insects, 
including honeybees, with audio signals. The re- 
searchers conducted an extensive evaluation of dif- 
ferent feature sets from audio analysis and machine 
learning algorithms for the audio classification task. 
Combining multiple different classifiers based on the 
same feature set was shown to be ineffective while en- 
sembles over different feature sets obtained from dif- 
ferent signal representations performed better. 

3 Hardware 

The EBM literature, briefly reviewed in the previous 
section, indicates that CV remains relatively unex- 
plored in EBM. However, our position is that CV is 
ready for inclusion in standard EBM sensor suites. 
CV is one of the virtual sensors of BeePi, a multi- 
sensor, solar-powered EBMD [5]. Each EBMD con- 
sists of a Raspberry Pi (RPi) computer, a camera, a 
solar panel, a temperature sensor, a battery, a hard- 
ware clock, a set of microphones, and a solar charge 
controller. 

The exact hardware components of BeePi, shown 
in figure (2), include an RPi 3 Model B with 1GB 
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RAM, an RPi T-Cobbler, a full-size breadboard for 
sensor integration, a waterproof DS18B20 digital tem- 
perature sensor, an RPi Camera Board (CB), and a 
USB microphone hub where three microphones are 
plugged in. The wired microphones are connected to 
the RPi and placed into small holes drilled in hive 
walls to collect audio samples from inside the hive. 
The holes are covered with mesh nets on the inside to 
prevent propolisation of microphones by the bees. 



Figure 2: BeePi hardware components: 1) battery; 2) 
RPi; 3) RPi camera board; 4) car charger; 5) bread- 
board; 6) solar charge controller; 7) solar panel wires 

Unlike the two previous hardware realizations of 
the BeePi EBMD described in [5, 6], the current hard- 
ware realization has the camera not only placed under 
a plastic cover but also attached to a wooden plank 
for improved balance and stability against wind. The 
plank is attached with screws and metallic brackets to 
the hive super where the BeePi hardware components 
reside, as shown in figure (3). The camera is further 
protected from the elements by a wooden box open at 
the bottom and attached to a migratory hive lid with 
screws and metallic brackets. When the lid is placed 
on the super with the BeePi hardware, as shown in 
figure (4), the box protects the camera against the 
elements from above and the four sides. 



Figure 3: Camera (1) under plastic cover attached 
to wooden plank (2); plank is attached with metal- 
lic brackets (3) to box with BeePi hardware; camera 
looks down on landing pad (4) 


Renogy 50 watts 12 Volts monocrystalline so- 
lar panels are used for harvesting solar energy into 
UPG 12V 12 Ah F2 sealed lead acid AGM deep-cycle 
rechargeable batteries via Renogy 10 Amp PWM solar 
charge controllers. It takes approximately 25 minutes 
to wire a BeePi EBMD for deployment. Figure (5) 
shows the first author wiring two BeePi units for field 
deployment in May 2016. 



Figure 4: Hive lid (1) with protection box (2) on hive; 
solar panel (3) on ground 


Data collection is done on the RPi with the soft- 
ware written in Python 2.7.9. The collected data are 
saved on a 32GB sdcard inserted into the RPi. When 
the system starts, three data collection threads are 
spawned. The first thread collects temperature read- 
ings every 5 minutes and saves them into a text file. 
The second thread collects 30-second wav recordings 
every 5 minutes. The third thread captures pictures of 
the hive’s landing pad every 5 minutes and saves them 
as PNG. The bee counting method presented in this 
article works with these images. A cronjob monitors 
the threads and restarts them when necessary. 



Figure 5: Wiring BeePi for field deployment 


4 Landing Pad Localization 

The proposed bee counting method consists of two 
stages: landing pad localization and bee counting on 
the localized bee pads. In the first stage, landing pads 
are localized and cropped from images. In the second 
stage, bees are counted on localized pads. This section 
describes the pad localization stage. Bee counting is 
discussed in Section 5. 

A pad image, as seen in the left image in figure (6), 
is converted to grayscale and the average brightness 
is computed. If the brightness is below a threshold, 
which in the current implementation is set to 190, the 
original image is converted from RGB to YCrCb, where 
Y is the luminance component and Cr and Cb are the 
blue difference and red difference chroma components, 
respectively. The converted image is split into the Y, 
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Figure 6: Left: original image; right: image with ad- 
justed brightness 

Cr, and Cb channels, and the histogram is equalized for 
each channel. The channels are merged into one image 
and the image is converted back to RGB, as shown in 
the right image in figure (6). 



Figure 7: An image of a skewed pad 

As shown in figure (7), some captured images have 
skewed landing pads. The skew may be caused by 
strong wind or by the physical manipulation of the 
beehive by a beekeeper. Therefore, the next step in 
the pad localization is to determine the skew angle of 
the pad and, if there is a skew, to rotate the image 
accordingly. Pad skew detection is based on the ob- 
servation that the grass around the landing pad has 
more edges than the landing pad itself even when the 
bees are present on it. One morphological operation 
that captures this texture difference is corner detec- 
tion, because the grass is naturally expected to have 
more corners than the landing pad. 



Figure 8: Left: corners in original image; right: true 
corner pixels 


Toward that end, the image with adjusted bright- 
ness is processed with the Harris Corner detector [15] 
and the maximum pixel value is recorded as v max (see 
the left image in figure (8)). A clone of the image 
with detected corners is eroded. The image with de- 
tected corners is compared with the image with eroded 
corners and only unmodified pixels are retained and 
thresholded at or above v max . Then AND operation is 


applied to the image with detected corners and the 
image with thresholded corners, to delete all modified 
pixels. In figure (8), the left image shows the image 
with detected corners and the right image shows the 
image after the AND operation. 

After the corners are detected, a hash map, 
M r (i) — >> (co, ci, ..., c n ), is created where i is a row 
number and Cj , 0 < j < n are the columns of the de- 
tected corners in row r. Another hash map, M c (i) — >• 
(ro, ri, ..., r n ), is created where i is a column number 
and r j , 0 < j < n are the rows of the detected corners 
in c. 

The maps M r and M c are used to construct ver- 
tical and horizontal frequency projections, i.e., two 
corner frequency histograms H r and H C1 respectively. 
In H ri the bin numbers correspond to rows and the 
values of the bins are the counts of detected corners 
in those rows; in iL c , the bin numbers correspond to 
columns and the values of the bins are the counts of 
true corners in those columns. The images in figures 
(9) and (11) show H r and i7 c , respectively, computed 
for the image in figure (7). 



100 150 200 250 300 


Figure 9: Row histogram H r of image in figure (7) 



Figure 10: Curve fit r{pc) over H r in figure (9) 

The histograms H r and H c are used to detect de- 
clines in detected corner counts. Since the landing 
pads have fewer corners than grass, these declines in- 
dicate the upper corner points of the landing pad. To- 
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Figure 11: Column histogram H c of image in figure 
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Figure 12: Curve fit c(x ) over H c in figure (11) 

ward that end, the histograms are converted into 2D 
scatter plots and the outliers in the plots are removed 
with the LI distance in which the distance between 
two points is the sum of the absolute differences of 
their Cartesian coordinates. A curve fitting procedure 
is applied to the corner counts to model the distribu- 
tion of corner pixels with a 4-th degree polynomial. 
Equation (1) gives the curve for rows; equation (2) 
gives the curve for columns. The coefficients of x 4 and 
x 3 are very small in both curves, and are discarded. 

r(x) = 0.0025x 2 + 0.017x + 24.7 (1) 

c(x) = -O.OOOlx 2 + 0.004a; + 6.4 (2) 

The curves in equations (1) and (2) are used to 
locate the landing pad in a given image. Row-based 
localization uses r(x) in equation (1) shown in fig- 
ure (10). A window W[ of length 10 is created 
and shifted by 5 units left to right along the val- 
ues of r(x). Let Mi be the median of W[ and let 
A \ — Mi — Mi _ ;|. If < 0, the curve’s slope is de- 
creasing, which indicates the start of the landing pad. 
Let I r = {i |A£ < 0} and let Qi r be the first quar- 
tile point of I r . The first quartile point approximates 
the upper ordinate of the landing pad. Since the pad 
is approximately fifty pixels wide, given the distance 


of the camera from the landing pad, the lower ordi- 
nate of the landing pad in the image is computed as 
Qi r + 50. Figure (13) shows the row-based cropping of 
the landing pad region obtained with this procedure 
applied to the image in figure (7). 



Figure 13: Row-based cropping of landing pad region 
in image in figure (7) 

After row-based pad localization, a similar pro- 
cedure is executed for column-based pad localization 
using c(x) in equation (2) shown in figure (12). From 
the curve one can see that there is a decline followed 
by a rise. A decline indicates a steady drop in the 
number of column corner counts, which, in turn, in- 
dicates the horizontal start of the landing pad. On 
the other hand, a rise indicates the horizontal end 
of the landing pad. A window W? of length 10 is 
created and shifted by 5 units left to right along the 
values of c(x). Let M? be the median of WJ. Let 
A f = M^ — Mf_ x . If A i < 0, the curve’s slope is de- 
creasing. Let I c = {i|A? < 0}. Let M £ be the median 
of I c . This value represents the left column where the 
landing pad is likely to start, i.e., the left horizontal 
start of the landing pad. Let J c = {j \Aj > 0} and let 
Mg be the median of J c . Then M e c is the right column 
index of the landing pad, i.e., the horizontal end of 
the landing pad. Figure (14) shows the column-based 
cropping applied to the image in figure (13). 


Figure 14: Column-based cropping of landing pad re- 
gion in image in figure (13) 

As figure (14) indicates, there may be some 
grass left in the image after column-based cropping. 
To eliminate these grass spots, another iteration of 
column-based localization described above is executed 
after column-based cropping. Specifically, the corner 
detection algorithm is applied to the cropped image, 
as shown in figure (15), and the curve fitting proce- 
dure is applied to the corner counts to model the dis- 
tribution of corner pixels with a 4-th degree polyno- 
mial, as shown in figure (16). 



Figure 15: Detected corners in image in figure (14) 

The left and right segments of the polynomial 
curve are ignored after the leftmost and rightmost 
non-zero frequency counts. For example, in figure 16, 
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Figure 16: Fourth degree polynomial fit over corner 
distribution after 2 nd iteration of column-based pad 
localization applied to image in figure 15 


the left falling segment can be safely ignored from 0 
to the last non-zero column frequency count. The 
polynomial rise on the right, on the other hand, indi- 
cates the presence of grass. To detect the exact end of 
the landing pad, a window of 10 points W{ is created 
and shifted left to right 5 units at a time along the 
curve values. Let Mi be the median of the window 
points and let A i = {i\Mi — < 6}, 6 = 0.105 

and A r = {i\Mi — > 6}, 6 = 0.105. If m is 

the middle of the image, two indices k and j are se- 
lected such that k = min{i\i G A/ A i < m} and 
j = max{i\i G A r A i > m}. The middle columns of 
Wk and Wj mark the left and right horizontal bound- 
aries of the landing pad, respectively. 



Figure 17: Reduced landing pad after 2 nd iteration 
of column-based plan localization applied to image in 
figure (14) 



Figure 18: Corners detected in right half of cropped 
image in figure (14) 

The final steps in landing pad localization are skew 
angle detection and image rotation. Pad skew detec- 
tion is similar to the method we previously developed 
to detect text skew in nutrition labels on grocery prod- 
ucts [16]. The input to the skew angle detection mod- 
ule is a horizontally and vertically cropped landing 
pad image like the one in Figure (17). The image is 
divided into the left and right halves. Another round 
of corner detection is executed on both halves. If the 
number of detected corners in either half is below a 
threshold the image is discarded. For example, the left 


half of the image in figure (14) does not have enough 
corners. However, the right half of the image in figure 
(14) does have some corners, as shown in figure (18). 



Figure 19: Skew angle detection in image in figure 

( 15 ) 

A list of rows, Ri, is computed in which there are 
no corners detected. Let min(Ri) be the minimum 
row with no corners detected with row 0 being the top 
row in the image. Let max(Ri) be the maximum row 
with no corners detected with row 0 being the bottom 
row in the image i.e scanning the image from bottom 
to top. If | min(Ri) — max(Ri)\ <6,6 = 3, the pad 
has no skew and does not need to be rotated. Oth- 
erwise, a column-based scan of the corner counts left 
to right indicates whether the corner counts increase 
left to right or decrease. A line is drawn between the 
points (min(Ri), 0) and ( max(Ri),w — 1 ), where w is 
the width of the image. This virtual line is depicted by 
a yellow dash line in figure (19). If the corner counts 
increase, the skew angle is Z(0,0), ( max(Ri),w — 1 ), 
(max(Ri), 0). If the counts decrease, the skew an- 
gle is Z(0, w — 1 ), (max(Ri), 0), ( max(Ri),w — 1 ). In 
figure (19), the hypotenuse [(0, 0), (max(Ri), w — 1 )] 
of Z(0, 0), ( max(Ri),w — 1 ), (max(Ri), 0), is depicted 
with a red dashed line. Figure 20 shows the rotated 
image. 



Figure 20: Pad image in figure (17) rotated counter- 
clockwise to eliminate skew 


5 Bee Counting 

After the landing pad is cropped from the image, the 
bees are counted on it. This section presents two bee 
counting algorithms we have developed so far. In Sec- 
tion 6, the performance of these algorithms is evalu- 
ated on a sample of 793 images, which corresponds 
to the daytime periods of three consecutive calendar 
days. 


5.1 Algorithm 1 

The first bee counting algorithm is based on the ID 
Haar Wavelet Transform (ID HWT) [17]. In the ID 
HWT, a signal is a vector in R n , n = 2 k , k G N. Fol- 
lowing the formalization in [18], let be a 2 k x 2 k 
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for computing k scales of the ID HWT. This matrix 
can be effectively computed from the n canonical base 
vectors of R n . If x = (xo, ...,x 2 k_ 1 ) is a signal in i? n , 
then y is the k th -scdle ID HWT of x is defined in 
equation (3). 


. X T = y (3) 


The values of the column vector y are defined in (4). 


K 


5 c 0 


O-i 
C 2* 


- 1 -!) 


( 4 ) 


In equation (4), = ju(y) and cj is the coefficient 

of the i th Haar wavelet at scale j. 




Figure 21: Up- down spikes 

HWTs are used to detect significant changes in sig- 
nal values [19]. We have previously proposed to clas- 
sify signal changes as wavelet spikes [20] . Specifically, 
we propose four types of spikes to detect and classify 
signal changes: up-down triangle, up-down trapezoid, 
down- up triangle, and down- up trapezoid. Figure (21) 
shows up-down triangle and trapezoid spikes; figure 
(22) shows down-up triangle and down-up trapezoid 
spikes. In both figures, the lower graphs represent 
possible values of the corresponding Haar wavelets at 
a given scale k. The difference between up-down and 
down-up spikes is the relative positions of the climb 
and decline segments. In trapezoid spikes, flat seg- 
ments are always in between the climb and decline 
segments. Formally, a spike S is a 9-tuple whose ele- 
ments are real numbers given in (5). 


S = (u s ,u e ,a, f s , f e ,~t,d, s ,d e , P) 


( 5 ) 


The first two elements of the 9-tuple, u s and u e , 
are the abscissae of the beginning and end of the 
spike’s climb segment, respectively, when the wavelet 
coefficients of the ID HWT increase. This happens 
only if the samples taken from the analyzed signals 
also increase. If cu ^ and cu ^ are the k th scale 
wavelet coefficient ordinates at u s and u e , respectively, 
then the sharpness of the spike’s climb is measured as 



Figure 22: Down-up spikes 

a = tan~ 1 {u e — u s + 1 ,cu^ — cu ^), which is the 
third element of the 9-tuple. 

The decline segment of the spike S in equation (5) 
is characterized by the seventh, eighth, and ninth el- 
ements of the 9-tuple: d s , d e , and /?, where d s and 
d e are the abscissae of the beginning and end of the 
spike’s decline segment, respectively, when the wavelet 
coefficients decrease. The decline occurs when the sig- 
nal’s samples decrease. If cd^ and cd^ are the k th 
scale wavelet coefficient ordinates at d s and d e , re- 
spectively, then the sharpness of the spike’s decline is 
measured as /3 = tan~ 1 (d e — d s + 1, cd ^ — cd ^). 

In trapezoid up-down or down- up spikes, the flat 
segment of the spike in equation (5) is character- 
ized by the fourth, fifth, and sixth elements of the 
9-tuple: / s , / e , and 7, where f s and f e are the ab- 
scissae of the beginning and end of the spike’s flat 
segment, respectively, over which the wavelet coeffi- 
cients either remain at the same ordinate or have mi- 
nor ordinate fluctuations. If cf ^ and cd^ are the 
k th scale wavelet coefficients corresponding to f s and 
/ e , respectively, the spike’s flatness is estimated as 
7 = tan~ 1 {f e — f s + l,c/i^ — c/i* 0 ). The absolute 
value of 7 is close to 0. 

Given an image, spikes are computed for each row. 
Of course, depending on the application domain, col- 
umn spikes can also be computed. When spikes are 
computed for row r, the column indices of the ac- 
tual pixels covered by each spike at scale j are com- 
puted by the formula in equation (6), where s and e 
are the positions of the starting and ending wavelet 
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coefficients in the ID HWT at scale j, respectively. 
For up-down spikes, s = u s and d = d e , whereas, for 
down- up spikes, s = d s and e = u e . 

p(j , s, e) = {z|2 J • s < i < 2 J • (e + 1) — 1} (6) 

Let n be the number of rows in an image and let 
U r be the set of up-down spikes in row r, where 0 < 
r < n — 1. The set of pixel columns in row r covered 
by the up-down spikes in r is given in equation (7), 
where j is a scale and s z and e z are the beginning 
and end positions of an up-down spike Z in row r, 
respectively. 

Z u,j = U PU> s z,e z ) (7) 

zeu r 

Let D r be the set of down-up spikes in row r, 0 < 
r < n — 1. The set of pixel columns in r covered by 
the down- up spikes is given in (8), where j is a scale 
and s z and e z are the beginning and end positions of 
a down-up spike Z in row r, respectively. 

z dj = IJ PU,s z ,e z ) (8) 

zeD r 

The number of unique column pixels covered by 
the up-down and down-up spikes in row r is given in 
equation (9). The formula in (10) gives the actual 
number of pixels covered by the up-down and down- 
up spikes in an image / with n rows. 

7 = (9) 

Xj{I ) = J2\Zj\ (10) 

r=0 

For example, consider three 16 x 16 images in fig- 
ure (23). Suppose it is required to separate bee pixels 
from non-bee pixels in the original bee image on the 
left. The bee pixels are those covered by up-down and 
down-up spikes in each row. The middle image in fig- 
ure (23) shows the only up-down spike detected in row 
8 after a single scale of the ID HWT. The right image 
in figure (23) shows the only down-up spike detected 
in row 8. Both spikes are triangle ones. 


Figure 23: Bee image (left); up-down triangle spike 
in row 8 (middle); down- up triangle spike in row 8 
(right); up-segments are red; down-segments are blue 

The up-down spike, S u d, shown in the middle im- 
age of figure (23), is defined in (11). Per notation 


in equation (5), the spike’s climb segment starts at 
u s = 3 and ends at u e = 4. The angle of the climb 
is a « 7r/4. Since this is a triangle spike, the flat 
segment is characterized as f s = f e = — 1 and 7 = 0. 
The spike’s decline segment starts at d s = 5 and ends 
at d e = 5. The decline’s angle is /5 ~ 7r/3. 

S u d = (3, 4, tt/4, -1, -1, 0, 5, 5, tt/3) (11) 

The down-up spike, Sdu, shown in the right image 
of figure (23), is defined in (12). The spike’s climb 
segment starts at u s = 3 and ends at u e = 4. The 
angle of the climb is a « 7r/4. Since this is also a 
triangle spike, the flat segment is characterized as f s = 
f e = — 1 and 7 = 0. The spike’s decline segment starts 
at d s = 1 and ends at d e = 2. The decline’s angle is 
P « 7r/3. Note that since this is a down-up spike, the 
abscissae of the start and end of the decline segments 
are to the left of the abscissae of the start and end of 
the climb segment. 

= (3, 4, tt/4, -1, -1, 0, 1, 2, tt/3) (12) 

Using equations (6), (7), and (8), the pixel 
columns of the up-down and down- up spikes in (11) 
and (12), respectively, are given in equations (13) and 
(14), respectively. 

z u,i = {p( 1>3, 5)} = {i|6 <i < 11} (13) 

7y = 0?(l,l,4)} = {i|2<i<9} (14) 

Using equation (9), the set of pixel columns cov- 
ered by the two spikes in row 8 in the left image of 
figure (23) is given in equation (15). 

From (13), (14), and (15) one can compute the 
actual bee pixels in row 8, i.e., {10, 11} U {2, 3, 4, 5} U 
{6,7, 8,9} = {i\2 < i < 11}. If the number of scales 
is j = 1, the number of bee pixels in the left image of 
figure (23) is given in equation (16). 

15 

X 1 (I)^\Z r 1 \=U (16) 

r = 0 

The pseudocode of the bee counting algorithm is 
given in figure (24). The algorithm takes an image / 
such as the image in figure (7), the normalizer TV, and 
the number of scales j of the ID HWT. The algorithm 
also takes the thresholds for the angles of up-down, 
flat, and down- up spikes, i.e., cq 7, /T In the current 
implementation, a = (3 = 60° and = 5°. 

5.2 Algorithm 2 

Algorithm 2 starts by image brightness adjustment 
and continues with the iterative reduction of the im- 
age to the actual landing pad to compensate for poten- 
tial inaccuracies in landing pad localization outlined 
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1. procedure countBeesl(/, N, j, a, 7 , P) 

2. L = localizeLandingPad(/); 

3. gaussianBlur(T); 

4. pyramidMeanShiftFilter(T); 

5. maxRGBFilter(T); 

6 . bleachBluePixels(T); 

7. convert ToGray scale (T); 

8 . return [Xj/N\; 


Figure 24: Bee counting algorithm 1 


in Section 4. As in algorithm 1 described in Section 
5.1, bee counts are computed by dividing the total 
number of bee pixels by the average number of pix- 
els occupied by individual bees obtained from camera 
calibration experiments. 

Image brightness is maximal when the sun is di- 
rectly above the beehive. When the sun is obscured by 
clouds, captured images become darker. To minimize 
the adverse effect of too much brightness or darkness, 
image brightness is dynamically adjusted to lie in in 
the range (45, 95), i.e., the brightness index should 
be greater than 45 but less than 95. This range was 
experimentally found to yield optimal results. 



Figure 25: Adjusting brightness in image with local- 
ized landing pad 

Figure (25) illustrates how brightness adjustment 
improves bee counts. The upper image on the right 
in figure (25) shows a green landing pad localized in 
the left image without adjusted brightness. The lower 
image on the right in in figure (25) shows a green pad 
identified in the same image with adjusted brightness. 
Only four bees were identified in the upper image on 
the right whereas in the lower image eight bees were 
identified, which is closer to the twelve bees found in 
the original image by a human evaluator. 


three steps of grass area reduction are shown in fig- 
ure 26. The first step is to convert the input from 
RGB to HSV according to equations (17) and (18). In 
equation (17), R ' , G', and B' are the R, G, B val- 
ues normalized by 255, Cmax = max{R\ G', B'} and 
A = Cmax ~ max{R' , G ' , B'}. The value of V is set 
to G maa ,. 

In the actual implementation, the format conver- 
sion is accomplished with the cvtColorO method of 
OpenCV. The inRangeO method of OpenCV is sub- 
sequently applied to identify the areas of green or 
white pixels, the two colors in which the landing pads 
of our beehives are painted. Noise is removed through 
a series of erosions and dilations. The white pixels in 
the output image represent green or white color in the 
actual image and the black pixels represent any color 
other than green or white. 



Figure 26: Landing pad identification steps: 1 ) HSV 
conversion; 2) color range identification; 3) noise re- 
moval 

To further remove noise from the image and to 
reduce it as closely as possible to the actual landing 
pad, contours are computed with the f indContours () 
method of OpenCV and a bounding rectangle is found 
for each contour. The bounding contour rectangles 
are sorted in increasing order by the Y coordinate, 
i.e., increasing rows. Figure (27) shows the bounding 
rectangles for the contours computed for the output 
image of step 3 in figure (26). 



Figure 27: Bounding rectangles of found contours 
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To compensate for inaccuracies in landing pad lo- 
calization described in Section 4, which may leave 
some small grass area above or on either side of the 
actual landing pad, algorithm 2 reduces the grass area 
present in the image of a localized landing pad. The 


Data analysis indicates that if the area of a con- 
tour is at least half the estimated area of the landing 
pad, the contour is likely to be part of the actual land- 
ing pad. On the other hand, if the area of a contour is 
less than 20 pixels, that contour is likely to be noise. 
In the current implementation of algorithm 2, the es- 
timated area of the green landing pad is set to 9000 
pixels and the estimated area of the white landing pad 
is set to 12000 pixels. These parameters can be ad- 
justed for distance, depending on how far the BeePi 
camera is from the landing pad. Using this pixel area 
size filter, the approximate location of the landing pad 
is computed by scanning through all the contours in 
the sorted list and finding the area of each contour. 
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If the area is at least half the estimated size of the 
landing pad of the appropriate color, the Y coordi- 
nate of the contour rectangle is taken to be the mean 

Y coordinate and the scanning process terminates. 

If the contour’s area is between 20 and half the 
estimated landing pad area, the Y coordinate of the 
contour is saved. Otherwise, the current contour is 
skipped and the next contour is processed. When the 
first contour scan terminates, the mean Y coordinate, 
/i(y), is calculated by dividing the sum of the saved 

Y coordinates by the number of the processed contour 
rectangles. 

After fi{Y) is computed, the second scan of the 
sorted list of contour rectangles is performed to find 
all contours whose height lies in [ji(Y) — H , fi(Y) + H], 
where H is half the estimated height of the landing 
pad for the appropriate color. In the current imple- 
mentation, H is set to 20 for green pad images and to 
25 for white pad images. The parameter H is hive de- 
pendent, because the exact camera location may differ 
slightly from hive to hive. For example, if the camera 
is placed closer to the landing pad, then H will have a 
higher value and if the camera is placed far from the 
landing pad, H will have a lower value. 

Bounding rectangles are computed after the sec- 
ond scan to enclose all points in the found contours, 
as shown in figure (27). To verify whether the im- 
age has been sufficiently reduced to the actual landing 
pad, the area of the bounding rectangle is computed. 
If the area of the bounding rectangle is greater than 
the estimated area of the landing pad, the bounding 
rectangle may still contain noise, in which case an- 
other scan is iteratively performed to remove noise by 
decreasing H by 2 to 4 units. In most of the cases, 
this extra scan is not needed, because the landing pad 
is accurately found after the first two scans. Figure 
(28) illustrates the three steps of the contour analysis 
to reduce the image to the actual landing pad. 





Figure 28: Contour analysis: 1) 1 st contour scan; 2) 
2 nd contour scan; 3) pad cropping 


pad image are arbitrarily set to 255 to facilitate the 
upcoming bee counts. Figure (29) shows the output 
of background separation. In figure (29), the green 
background is converted to white and the foreground 
to black. Since noise may be introduced, the image 
is de-noised through a series of erosions and dilation 
with a 2 x 2 structuring element. 


Figure 29: Background and foreground separation 

To count bees in the image after background sep- 
aration, the image is converted to grayscale and the 
contours are computed again. Data analysis suggests 
that the area of an individual bee or a group of bees 
vary from 20 to 3000 pixels. Therefore, if the area of 
a contour is less than 20 pixels or greater than 3000 
pixels, the contour is removed. The area of one in- 
dividual bee is between 35 and 100 pixels, depending 
on the distance of the camera from the landing pad. 
The green landing pad images were captured by a pi 
camera placed approximately 1.5 m above the landing 
pad with the average area of the bee being 40 pixels. 
On the other hand, the white landing pad images were 
captured by a pi camera placed approximately 70 cm 
above the landing pad where the average area of an 
individual bee is 100 pixels. To find the number of 
bees in green landing pad images, the number of the 
foreground pixels, i.e., the foreground area, is divided 
by 40 (i.e., the average bee pixel area on green land- 
ing pads), whereas, for the white landing pad images, 
the foreground area is divided by 100 (i.e., the average 
bee pixel area on white landing pads). The result is 
the count of bees in the image. In the upper image 
in figure (30), five bees are counted by the algorithm. 
The lower image in figure (30) shows the found bees 
in the original image. 



After the pad is cropped, foreground and back- 
ground pixels are separated on the basis of color. 
Specifically, in the current implementation of algo- 
rithm 2, for green landing pads, the background is 
assumed to be green and the foreground, i.e., the bees, 
is assumed to be yellow; for white landing pads, the 
background is assumed to be white and the foreground 
is assumed to be yellow. All pixels with shades of 
green or white are set to 255 and the remaining pixels 
are set to 0. Three rows of border pixels of the landing 


Figure 30: Omnidirectional bee counting 

The pseudocode of algorithm 2 is given in figure 
(31). The first input parameter to this algorithm is 
the image I like the image in figure (20). The param- 
eters Ci and C u specify the lower and upper pixel ar- 
eas for the contours. In the current implementation, 
Ci = 20 and C u = 3000. The last parameter A^ ee 
specifies the average pixel area of an individual bee. 
For green pads, A^ ee = 40; for white pads, Ab ee = 100. 
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1. procedure countBees2 (/, C/, C u , Ab ee ) 

2. A = adjust Brightness (I); 

3. H = convert ToHSV (A) ; 

4. P = cropLandingPad(H); 

5. B = separateBackground(P); 

6. CONS = identifyContours(B); 

7. AREA = 0; 

8. for each c in CONS 

9. if (area(c) Ci and area(c) Cu) 

10. AREA += area(c); 

11. end if 

12. end for 

13. BeeCount = [AREA/Ab ee \; 

14. return BeeCount; 


Figure 31: Bee counting algorithm 2 

6 Evaluation 

To compare the accuracy of the two algorithms, 793 
images of landing pads with bees on them were se- 
lected. Three human observers counted the bees on 
the pad in each image. The human counts served 
as ground truth. Both bee counting algorithms were 
evaluated on this sample of images. Table (1) sum- 
marizes the evaluation results. The first column gives 
the names of the two algorithms used. The numbers 
in the columns are error margin values that give the 
allowed margin of error between the human and com- 
puter counts for each image. For example, the value of 
the second column is equal to 3. Thus, the value 63.18 
is the percent of images where bee counting algorithm 
2 approached the human counts of the same images 
within a margin of 3, i.e., the absolute difference be- 
tween the bees counted by the algorithm and the hu- 
man count did not exceed 3. Similarly, the two values 
in the last column indicate that the first bee count- 
ing algorithm approached the human counts within 
a margin of 15 on 82.10% of the images, whereas the 
second bee counting algorithm approached the human 
counts within a margin of 15 on 90.29% of the images. 


Algo/Error 

3 

5 

7 

10 

15 

Algo 1 

39.97 

53.72 

63.56 

73.77 

82.10 

Algo 2 

63.18 

73.77 

80.10 

84.99 

90.29 


Table 1: Evaluation of two bee counting algorithms: 
column names are error margin values; other values 
are percentages 


V c 4- * M* 


Figure 32: Inaccurately localized landing pad 

The results in table (1) indicate that the second al- 
gorithm outperformed the first bee counting algorithm 


for each error margin. This relative underperformance 
can be attributed to the fact that, unlike algorithm 2, 
algorithm 1 does not compensate for inaccurate land- 
ing pad localization by reducing the size of the image 
to the actual landing pad. Figure (32) shows a typical 
image on which algorithm 2 outperforms algorithm 1. 
Note that the upper part of the image contains a siz- 
able area of grass that algorithm 2 eliminates in steps 
4 and 5 of the pseudocode in figure (31). Algorithm 
1, on the other hand, detects quite a few spikes in the 
grassy area above the actual landing pad, which re- 
sults in many false positives. When the landing pad is 
localized accurately both algorithms perform equally 
well. False positives were also caused by occasional 
shadows, leaves, or blades of grass wrongly counted 
as bees. 

7 Conclusion 

The investigation reported in this article continues 
our research on CV algorithms for omnidirectional 
bee counting [5, 6, 20]. The method for omnidi- 
rectional beecounting improves the previously pro- 
posed methods in three respects. First, the presented 
method completely automates landing pad localiza- 
tion whereas in our previously proposed methods the 
approximate landing pad coordinates were hardcoded 
into our program. Second, the presented method was 
implemented and evaluated in situ on a hardware 
device assembled with off-the-shelf hardware compo- 
nents. In our previous work, only image capture ran in 
situ whereas the implemented bee counting methods 
were evaluated in the laboratory on a standard desk- 
top computer [5]. Finally, bee counts are improved 
through automatic detection and elimination of land- 
ing pad skews, which had a negative impact on bee 
counts in our previous experiments [6]. The obtained 
results indicate that CV can be included in standard 
EBM sensor suites. We hope that CV will help re- 
searchers, practitioners, and citizen scientists to build 
more reliable EBMDs and to promote the grassroots 
development of regional, national, and international 
cyberinfrastructures for capturing, sharing, and ana- 
lyzing EBM data. 
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